%   brk2headmsk.m - outputs the mri scalp mask 0's outside the head and 1's
%   inside the head
%   Chris Whalen
%   University of Illinois
%   12/15/05 - Date Created
%   12/15/05 - Last Modified


% N.B. - later this can be integrated with brk2msk to save computational
% time if need be.

function [vol_mask] = brk2headmsk(file,mri_dim,vol,thresh)


fprintf('\n\n     brk2headmsk \n');
  
if isempty(thresh) % quick fix for not having threshold
    thresh = 0.05;
end

% Reading and Initializing vol=0
if vol == 0
    fid = fopen(file.brik);
    data=fread(fid,'uint16');
    fclose(fid);
    data=data/(max(data));  % rescale intensities 0 to 1
    vol = reshape(data,[mri_dim.x mri_dim.y mri_dim.z]);    % create BRIK
    vol = filter3d(vol); % Spatial Filters to remove graininess 
    vol = smooth3(vol); % gaussian filtering with kernel 3 and sd 0.65 . Must be done here or else xyz_mask won't be smoothed
    clear data
end

% Initializing
vol_mask = zeros(mri_dim.x, mri_dim.y, mri_dim.z); 





% File Being Processed
fprintf('\nEstimated Segmentation Time: 0.2 minutes\n');
tic;
  

% Generating Mask Left/Right
fprintf('Generating L-R Mask ...  \n');
%fprintf('Generating L-R Mask ... \n  Estimated Time: 0.2 min \n');
%tic;
change = 0;                            
for k=1:mri_dim.z                                       % Z Coords
    for j=1:mri_dim.y                                   % Y Coords
        for i=1:mri_dim.x                         % X Right Bound Scalp Coords (must do at least slightly more than half way in case head is asymetrical)
            if  (vol(i,j,k) >= thresh) & (change == 0) %(sum(vol(i,j,k)) == lin_clust_size)
                 change = 1;
                 vol_mask(i:mri_dim.x,j,k) = 1;
            end               
        end
        change = 0;
        for i=mri_dim.x:-1:(1)             % X Left Bound Scalp Coords 
            if  (vol(i,j,k) >= thresh) & (change == 0) %(sum(vol(i,j,k)) == lin_clust_size)
                 change = 1;
                 vol_mask(i:-1:1,j,k) = vol_mask(i:-1:1,j,k) + 1;
            end               
        end
        change = 0;    
    end
end
%t_left_right = toc;
%fprintf('  Actual Time: %-2.1f min\n', t_left_right/60);
    

% Generating Mask Anterior/Posterior
fprintf('Generating A-P Mask ... \n');
%fprintf('Generating A-P Mask ... \n  Estimated Time: 0.2 min \n');
%tic;
change = 0;                                 
for k=1:mri_dim.z                                       % Z Coords
    for i=1:mri_dim.x                                   % Y Coords
        for j=1:(mri_dim.y)                         % Anterior Bound Scalp Coords 
            if  (vol(i,j,k) >= thresh) & (change == 0) %(sum(vol(i,j,k)) == lin_clust_size)
                 change = 1;
                 vol_mask(i,j:mri_dim.y,k) = vol_mask(i,j:mri_dim.y,k) + 1;
            end               
        end
        change = 0;
        for j=(mri_dim.y):-1:1             % Posterior Bound Scalp Coords 
            if  (vol(i,j,k) >= thresh) & (change == 0) %(sum(vol(i,j,k)) == lin_clust_size)
                 change = 1;
                 vol_mask(i,j:-1:1,k) = vol_mask(i,j:-1:1,k) + 1;
            end               
        end
        change = 0;    
    end
end
%t_ant_post = toc;
%fprintf('  Actual Time: %-2.1f min\n', t_ant_post/60);

% Generating Mask Superior/Inferior
fprintf('Generating S-I Mask ... \n');
%fprintf('Generating S-I Mask ... \n  Estimated Time: 0.2 min \n');
%tic;
change = 0;                           
for i=1:mri_dim.x                                       % Z Coords
    for j=1:mri_dim.y                                   % Y Coords
        for k=(mri_dim.z):-1:1             % Superior Bound Scalp Coords 
            if  (vol(i,j,k) >= thresh) & (change == 0) %(sum(vol(i,j,k)) == lin_clust_size)
                 change = 1;
                 vol_mask(i,j,k:-1:1) = vol_mask(i,j,k:-1:1) + 1;
            end               
        end
        change = 0;    
    end
end
%t_sup_inf = toc;
%fprintf('  Actual Time: %-2.1f min\n', t_sup_inf/60);

thresh_crit = 3; % criterion for what vol number must be before it should be included in mask (1-5). needs to be at least 2. 5 means it appears in every direction.
vol_mask = vol_mask./thresh_crit; % scaling to desired thresh_crit;
vol_mask = floor(vol_mask); % this places into 0's and 1's based on thresh_crit;


% xyz_mask = xyz_mask + 1; % Adding 1 to have both BRIK and ELP be in same dimensions
t_brk2msk = toc;
fprintf('  Actual Time: %-3.0f min\n', t_brk2msk/60);

fprintf('\n     brk2headmsk Done \n');



